Covid19, Japan

日本の新型コロナウィルスに関する情報が抱える大きな問題は、ソースごとに値が微妙にことなっていることです。原因は明白で、総元締めであるはずの厚生労働省が各自治からの患者データを集計せずに各自治体のプレスリリースやウェブサイトの数値を積上げたものを発表していることにあります。
本来であれば医療機関->市町村->都道府県->厚生労働省の順で患者のデータを全て吸い上げて集計すべきなのです。

Import

厚生労働省

厚生労働省のオープンデータのページから取得した「陽性確定者数」のデータ

df_mhlw <- "https://www.mhlw.go.jp/content/pcr_positive_daily.csv" %>% 
  readr::read_csv() %>% 
  dplyr::mutate(date = lubridate::as_date(`日付`),
                n = as.integer(`PCR 検査陽性者数(単日)`),
                cum = cumsum(n)) %>% 
  dplyr::select(date, n, cum)

df_mhlw

　

Covid19Japan

Covid19 Japanが独自に収集している陽性確定者単位のデータ。ソースとデータは全てGitHubにて公開されているが、データはJSON形式である点に注意。発表後に修正されたレコード（インスタンス）は削除されれずにステータスなどが変更されているだけなので、「レコード数 \(\neq\) 累計陽性確定者」である点に注意。

df %>% 
  dplyr::group_by(patientId) %>% 
  dplyr::summarise(n = n())

df %>% 
  dplyr::filter(patientId == "-1")

df %>% 
  dplyr::filter(patientStatus == "Deceased")

JAG Japan

ジャッグジャパン株式会社が独自に収集している陽性確定者単位のデータ。データはCSV形式にて公開されているが、GIS処理用のデータが含まれている点に注意。「レコード数 \(=\) 累計陽性確定者数」になっている。

df_jag <- "https://dl.dropboxusercontent.com/s/6mztoeb6xf78g5w/COVID-19.csv" %>% 
  readr::read_csv()

## Warning: Missing column names filled in: 'X52' [52], 'X53' [53], 'X54' [54]

## Warning: 610 parsing failures.
##   row              col expected              actual                                                               file
##  4484 市町村内症例番号 a double 神戸59/西宮36       'https://dl.dropboxusercontent.com/s/6mztoeb6xf78g5w/COVID-19.csv'
##  4485 市町村内症例番号 a double 神戸60/西宮37       'https://dl.dropboxusercontent.com/s/6mztoeb6xf78g5w/COVID-19.csv'
##  6054 市町村内症例番号 a double 川口40/さいたま50   'https://dl.dropboxusercontent.com/s/6mztoeb6xf78g5w/COVID-19.csv'
## 12486 市町村内症例番号 a double 神戸237/西宮63      'https://dl.dropboxusercontent.com/s/6mztoeb6xf78g5w/COVID-19.csv'
## 14454 市町村内症例番号 a double (271)のちに番号なし 'https://dl.dropboxusercontent.com/s/6mztoeb6xf78g5w/COVID-19.csv'
## ..... ................ ........ ................... ..................................................................
## See problems(...) for more details.

df_jag

　

Area Data

都道府県と地方区分に関するデータをGistで公開していますので、使いたい方はご自由に。

　

Data Wrangling

ソースによる違い

収集組織	累計感染者数	累計死者数	日付
厚生労働省	96948		2020-10-26
Covid19 Japan	98146	1726	2020-10-27
JAG Japan	93661		2020-10-26
NHK（参考）	98262	1733	2020-10-27

https://www3.nhk.or.jp/news/special/coronavirus/data-all/

Summarize

オリジナルのデータがどのようになっているかskimrパッケージを用いてサマライズしてみます。元がJSON形式なので、読み込んだ直後は殆どの変量（フィーチャー）が文字型になっているのと欠損値が多いことが分かります。

df %>% 
  skimr::skim()

Data summary

Name	Piped data
Number of rows	99726
Number of columns	23
_______________________
Column type frequency:
character	19
logical	3
numeric	1
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
patientId	0	1.00	1	8	0	98147	0
dateAnnounced	0	1.00	10	10	0	273	0
gender	12662	0.87	1	1	0	2	0
detectedPrefecture	0	1.00	3	15	0	49	0
patientStatus	95831	0.04	8	23	0	8	0
notes	51724	0.48	1	270	0	45257	1
mhlwPatientNumber	99277	0.00	1	11	0	434	0
prefecturePatientNumber	10702	0.89	5	20	0	89015	0
prefectureSourceURL	68470	0.31	5	224	0	3433	0
residence	20455	0.79	1	38	0	1421	0
sourceURL	637	0.99	1	239	0	7736	0
relatedPatients	89684	0.10	2	259	0	6167	0
knownCluster	97269	0.02	3	88	0	228	0
detectedCityTown	74545	0.25	2	22	0	663	0
cityPrefectureNumber	74817	0.25	1	34	0	24900	2
citySourceURL	87930	0.12	9	317	0	3623	0
deceasedDate	97990	0.02	10	10	0	223	0
deceasedReportedDate	98609	0.01	10	62	0	183	0
deathSourceURL	98706	0.01	14	123	0	624	0

Variable type: logical

skim_variable	n_missing	complete_rate	mean	count
confirmedPatient	0	1	0.98	TRU: 98146, FAL: 1580
charterFlightPassenger	99712	0	1.00	TRU: 14
cruisePassengerDisembarked	99715	0	1.00	TRU: 11

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
ageBracket	0	1	32.98	23.44	-1	20	30	50	100	▃▇▅▂▁

　

Transform

各変量（フィーチャー）を適切な形式に変換し、地域区分でも分析できるように都道府県データと結合します。

x <- df %>% 
  dplyr::mutate(dateAnnounced = lubridate::as_date(dateAnnounced),
                ageBracket = forcats::as_factor(ageBracket),
                gender = forcats::as_factor(gender),
                patientStatus = forcats::as_factor(patientStatus),
                residence = forcats::as_factor(residence),
                cluster = dplyr::if_else(!is.na(knownCluster), TRUE, FALSE)) %>% 
  dplyr::select(dateAnnounced, ageBracket, gender, detectedPrefecture,
                patientStatus, knownCluster, cluster) %>% 
  dplyr::left_join(prefs, by = c("detectedPrefecture" = "pref"))

x

x %>% 
  skimr::skim()

Data summary

Name	Piped data
Number of rows	99726
Number of columns	13
_______________________
Column type frequency:
character	2
Date	1
factor	9
logical	1
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
detectedPrefecture	0	1.00	3	15	0	49	0
knownCluster	97269	0.02	3	88	0	228	0

Variable type: Date

skim_variable	n_missing	complete_rate	min	max	median	n_unique
dateAnnounced	0	1	2020-01-15	2020-10-27	2020-08-10	273

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
ageBracket	0	1.00	FALSE	13	20: 23833, 30: 14960, -1: 12773, 40: 12330
gender	12662	0.87	FALSE	2	M: 48990, F: 38074
patientStatus	95831	0.04	FALSE	8	Dec: 1726, Hos: 1266, Hom: 316, Dis: 283
pcode	1159	0.99	FALSE	47	13: 30790, 27: 12272, 14: 8559, 23: 6031
都道府県	1159	0.99	FALSE	47	東京都: 30790, 大阪府: 12272, 神奈川: 8559, 愛知県: 6031
八地方区分	1159	0.99	FALSE	8	関東地: 52254, 近畿地: 19525, 九州地: 10905, 中部地: 9643
広域圏	7448	0.93	FALSE	8	首都圏: 52471, 近畿圏: 18965, 中部圏: 8228, 九州圏: 7592
通俗的区分	1159	0.99	FALSE	11	関東: 52254, 関西: 18965, 東海: 7882, 九州: 7592
fct_pref	1159	0.99	FALSE	47	Tok: 30790, Osa: 12272, Kan: 8559, Aic: 6031

Variable type: logical

skim_variable	n_missing	complete_rate	mean	count
cluster	0	1	0.02	FAL: 97269, TRU: 2457

　
文字型を因子型に変換するだけでも、変量（フィーチャー）内の水準の比率が大まかにつかめるようになります。
例えば、年代別の陽性判定者数は20代が最も多く、続いて、30代、40代と働き盛りの世代に多いことが分かります。また、都道府県別では、東京、大阪、神奈川、愛知と成っていますが、地方区分別では、以外にも関東、大阪の次に九州がきていることが分かります。

　

Visualizing

厚生労働省

df_mhlw %>% 
  ggplot2::ggplot(ggplot2::aes(x = date)) + 
    ggplot2::
                  ma = zoo::rollmeanr(n, k = 7L, na.pad = TRUE)

全体傾向

sec_scale <- 50

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced) %>% 
  dplyr::summarise(n = n()) %>% 
  tidyr::complete(dateAnnounced = seq.Date(from = min(dateAnnounced),
                                           to = max(dateAnnounced),
                                           by = "day"),
                  fill = list(n = 0L)) %>% 
  dplyr::mutate(cum = cumsum(n),
                ma = zoo::rollmeanr(n, k = 7L, na.pad = TRUE)) %>%
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced)) + 
    ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity",
                      fill = "dark green", alpha = 0.5, width = 1.0) + 
    ggplot2::geom_line(ggplot2::aes(y = cum / sec_scale),
                       colour = "dark green") +
    ggplot2::geom_line(ggplot2::aes(y = ma), colour = "dark red") + 
    ggplot2::scale_y_continuous(
      name = "陽性確定数（日別）・7日間移動平均（濃赤折線）",
      sec.axis = ggplot2::sec_axis(~ . * sec_scale,
                                    name = "累積陽性確定者数（折線）")
    ) +
    ggplot2::labs(title = paste0("Covid19, Japan（@", Sys.time(), "）"),
                  subtitle = "", x = "日付")

## Warning: Removed 6 row(s) containing missing values (geom_path).

plotly::ggplotly()

地方区分別

都道府県別では区分が多すぎるので八地方区分別の積上げ棒グラフを描いてみます。サマリで見たように九州地方が以外にも多いことが読み取れます。

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced, `八地方区分`) %>% 
  dplyr::summarise(n = n()) %>% 
  tidyr::complete(dateAnnounced = seq.Date(from = min(dateAnnounced),
                                           to = max(dateAnnounced),
                                           by = "day"),
                  fill = list(n = 0L)) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, fill = `八地方区分`)) + 
    ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity", alpha = 1.0) +
    ggplot2::scale_fill_brewer(palette = "Set1")

折線グラフで表示してみます。

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced, `八地方区分`) %>% 
  dplyr::summarise(n = n()) %>% 
  tidyr::complete(dateAnnounced = seq.Date(from = min(dateAnnounced),
                                           to = max(dateAnnounced),
                                           by = "day"),
                  fill = list(n = 0L)) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, colour = `八地方区分`)) + 
    ggplot2::geom_line(ggplot2::aes(y = n)) +
    ggplot2::geom_point(ggplot2::aes(y = n), size = 1) + 
    ggplot2::scale_colour_brewer(palette = "Set1")

　
さらに分かりやすくするために関東、中部、近畿、九州の四地区を除く他の地区をその他としてまとめてみます。

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::mutate(`八地方区分` = forcats::fct_collapse(`八地方区分`,
                                                `その他地方` = c("北海道地方",
                                                          "東北地方",
                                                          "中国地方",
                                                          "四国地方"))) %>% 
  dplyr::group_by(dateAnnounced, `八地方区分`) %>% 
  dplyr::summarise(n = n()) %>% 
  tidyr::complete(dateAnnounced = seq.Date(from = min(dateAnnounced),
                                           to = max(dateAnnounced),
                                           by = "day"),
                  fill = list(n = 0L)) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, fill = `八地方区分`)) + 
    ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity", alpha = 1.0) +
    ggplot2::scale_fill_brewer(palette = "Set1")

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::mutate(`八地方区分` = forcats::fct_collapse(`八地方区分`,
                                                `その他地方` = c("北海道地方",
                                                          "東北地方",
                                                          "中国地方",
                                                          "四国地方"))) %>% 
  dplyr::group_by(dateAnnounced, `八地方区分`) %>% 
  dplyr::summarise(n = n()) %>% 
  tidyr::complete(dateAnnounced = seq.Date(from = min(dateAnnounced),
                                           to = max(dateAnnounced),
                                           by = "day"),
                  fill = list(n = 0L)) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, colour = `八地方区分`)) + 
    ggplot2::geom_line(ggplot2::aes(y = n)) +
    ggplot2::scale_colour_brewer(palette = "Set1") +
    ggplot2::facet_wrap(~ `八地方区分`, ncol = 2)

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced, `八地方区分`) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::ungroup() %>% 
  tidyr::pivot_wider(names_from = `八地方区分`, values_from = n) %>% 
  dplyr::mutate_if(is.integer, list(~tidyr::replace_na(., 0L))) %>% 
  tidyr::pivot_longer(cols = -dateAnnounced,
                      names_to = "八地方区分", values_to = "n") %>% 
  dplyr::mutate(`八地方区分` = forcats::fct_relevel(`八地方区分`,
                                              "北海道地方", "東北地方",
                                              "関東地方", "中部地方",
                                              "近畿地方", "中国地方",
                                              "四国地方", "九州地方")) %>%
  dplyr::group_by(`八地方区分`) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, y = cum,
                               fill = `八地方区分`)) +
    ggplot2::geom_area(stat = "identity", alpha = 1,
                       position = ggplot2::position_stack(reverse = FALSE)) +
    ggplot2::scale_fill_brewer(palette = "Set3") +
    ggplot2::labs(title = paste0("Covid19（@", Sys.time(), "）"),
                  y = "累計人数")

　
九州が8月頃から急増しているのは、県別に見ると福岡と沖縄での急増が原因と分かります。

x %>% 
  dplyr::filter(`八地方区分` == "九州地方") %>% 
  dplyr::group_by(dateAnnounced, `都道府県`) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::ungroup() %>% 
  tidyr::pivot_wider(names_from = `都道府県`, values_from = n) %>% 
  dplyr::mutate_if(is.integer, list(~tidyr::replace_na(., 0L))) %>% 
  tidyr::pivot_longer(cols = -dateAnnounced,
                      names_to = "都道府県", values_to = "n") %>% 
  dplyr::mutate(`都道府県` = forcats::fct_relevel(`都道府県`,
                                              "福岡県", "佐賀県", "長崎県",
                                              "熊本県", "大分県", "宮崎県",
                                              "鹿児島県", "沖縄県")) %>%
  dplyr::group_by(`都道府県`) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, y = cum,
                               fill = `都道府県`)) +
    ggplot2::geom_area(stat = "identity", alpha = 1,
                       position = ggplot2::position_stack(reverse = FALSE)) +
    ggplot2::scale_fill_brewer(palette = "Set3") +
    ggplot2::labs(title = paste0("Covid19, 九州地方（@", Sys.time(), "）"),
                  y = "累計人数")

　

クラスタ比率

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(`八地方区分`, cluster) %>% 
  dplyr::summarise(n = n()) %>% 
  tidyr::pivot_wider(names_from = cluster, values_from = n) %>% 
  dplyr::mutate(ratio = (`TRUE` / (`TRUE` + `FALSE`) * 100) %>% round(1)) %>% 
  dplyr::rename(`地方` = `八地方区分`,
                `非クラスタ感染者` = `FALSE`, `クラスタ感染者` = `TRUE`,
                `クラスタ比率[%]` = ratio)

Modeling

時系列（TS）分析

日毎の陽性判定者数を時系列データ形式に変換します。判定者数がゼロの日は報告されていないので、最初の陽性判定者が出た日から日日次のカレンダーデータを作成して結合してます。なお、時系列データ形式の周期については日次データなので7日間を周期として設定しておきます（考え方によっては約1ヶ月ごとで4週間=28日を周期にするのもありかと思います）。

tmp <- x %>%
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced) %>% 
  dplyr::summarise(n = n()) %>% print() %>% 
  dplyr::rename(date = dateAnnounced) %>% 
  tidyr::complete(date = seq.Date(min(date), max(date), by = "day"),
                  fill = list(n = 0L)) %>% print() %>% 
  # tidyr::replace_na(list(n = 0L)) %>% print() %>% 
  dplyr::select(-date)

## # A tibble: 270 x 2
##    dateAnnounced     n
##    <date>        <int>
##  1 2020-01-15        1
##  2 2020-01-24        1
##  3 2020-01-25        1
##  4 2020-01-26        1
##  5 2020-01-28        3
##  6 2020-01-29        1
##  7 2020-01-30        4
##  8 2020-01-31        1
##  9 2020-02-01        2
## 10 2020-02-04        1
## # … with 260 more rows
## # A tibble: 287 x 2
##    date           n
##    <date>     <int>
##  1 2020-01-15     1
##  2 2020-01-16     0
##  3 2020-01-17     0
##  4 2020-01-18     0
##  5 2020-01-19     0
##  6 2020-01-20     0
##  7 2020-01-21     0
##  8 2020-01-22     0
##  9 2020-01-23     0
## 10 2020-01-24     1
## # … with 277 more rows

ts_x <- tmp %>%
  ts(frequency = 1)

tsw_x <- tmp %>%
  ts(frequency = 7)

tsm_x <- tmp %>%
  ts(frequency = 28)

# ts_x <- x %>% 
#   dplyr::filter(!is.na(fct_pref)) %>% 
#   dplyr::group_by(dateAnnounced) %>% 
#   dplyr::summarise(n = n()) %>% 
#   dplyr::left_join(date, ., by = c("date" = "dateAnnounced")) %>% 
#   dplyr::mutate(n = tidyr::replace_na(n, 0L)) %>% 
#   dplyr::select(-date) %>%
#   ts(frequency = 1)
#   
# tsw_x <- x %>% 
#   dplyr::filter(!is.na(fct_pref)) %>% 
#   dplyr::group_by(dateAnnounced) %>% 
#   dplyr::summarise(n = n()) %>% 
#   dplyr::left_join(date, ., by = c("date" = "dateAnnounced")) %>% 
#   dplyr::mutate(n = tidyr::replace_na(n, 0L)) %>% 
#   dplyr::select(-date) %>%
#   ts(frequency = 7)
# 
# tsm_x <- x %>% 
#   dplyr::filter(!is.na(fct_pref)) %>% 
#   dplyr::group_by(dateAnnounced) %>% 
#   dplyr::summarise(n = n()) %>% 
#   dplyr::left_join(date, ., by = c("date" = "dateAnnounced")) %>% 
#   dplyr::mutate(n = tidyr::replace_na(n, 0L)) %>% 
#   dplyr::select(-date) %>%
#   ts(frequency = 28)

時系列データに変換したものをプロットすると可視化の項でプロットした棒グラフと同じ形のグラフになることが分かります。

ts_x %>% 
  plot()

tsw_x %>% 
  plot()

tsm_x %>% 
  plot()

上記からトレンド（長期的傾向）を除いたグラフ。デフォルト指定なのでlag = 1

ts_x %>% 
  base::diff() %>% 
  plot()

tsw_x %>% 
  base::diff() %>% 
  plot()

tsm_x %>% 
  base::diff() %>% 
  plot()

トレンド、季節変動（周期変動）、非周期変動に分解した場合。frequency = 1では分解できない点に注意。

tsw_x %>% 
  stats::decompose() %>% 
  plot()

tsm_x %>% 
  stats::decompose() %>% 
  plot()

トレンドだけを抜き出してみる。

tsw_x %>% 
  stats::decompose() %>% 
  .$x %>% 
  plot(ylim = c(0, 1500))

par(new = TRUE)

tsw_x %>% 
  stats::decompose() %>% 
  .$trend %>% 
  plot(ylim = c(0, 1500), col = "dark green", lwd = 3)

par(new = FALSE)

tsm_x %>% 
  stats::decompose() %>% 
  .$x %>% 
  plot(ylim = c(0, 1500))

par(new = TRUE)

tsm_x %>% 
  stats::decompose() %>% 
  .$trend %>% 
  plot(ylim = c(0, 1500), col = "blue", lwd = 3)

par(new = FALSE)

# tsm_x %>% 
#   stats::decompose() %>% 
#   .$trend %>% 
#   plot()

スペクトル密度の推定。何を意味するのかよく分からない。

ts_x %>% 
  spec.pgram()

tsw_x %>% 
  spec.pgram()

tsm_x %>% 
  spec.pgram()

自己回帰によるスペクトラム解析

spectrum(ts_x, method="ar")

spectrum(tsw_x, method="ar")

spectrum(tsm_x, method="ar")

Infer

ARIMA Model

ARIMA（Auto Regressive Integrated Moving Average, 自己回帰和分移動平均）モデルによる予測を行ってみます。予測に必要なパラメータはステップワイズにより自動的に最適なものが選択されます。

ts_x %>% 
  forecast::auto.arima() %>% 
  forecast::forecast() %>% 
  plot()

tsw_x %>% 
  forecast::auto.arima() %>% 
  forecast::forecast() %>% 
  plot()

tsm_x %>% 
  forecast::auto.arima() %>% 
  forecast::forecast() %>% 
  plot()